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ABSTRACT 

In  the  numerical  integration  of  the  primitive  equations  in 
spherical  coordinates  using  explicit  time  finite  difference  schemes, 
the  maximum  stable  time  step  decreases  toward  the  poles  as  increasingly 
smaller-scale,  higher-frequency  inertia-gravity  waves  are  resolved. 
Arakawa  has  developed  two  methods,  a "three-point"  method  and  a 
"Fourier"  method,  to  modify  the  frequencies  of  these  waves  such  that, 
without  Fourier  filtering  or  reducing  the  longitudinal  resolution, 
they  are  stable  for  a time  step  based  on  a prescribed  constant  length 
scale.  Experiments  have  shown  that  the  three-point  method  is  sub- 
stantially faster  per  time  step  than  the  Fourier  method  but  requires 
a time  step  approximately  half  that  of  the  Fourier  method  for 
stability. 

In  this  paper  the  three-point  method  is  analyzed  and  compared  to  j 
the  Fourier  method.  It  is  shown  that  the  original  definition  of 
the  latitudinal ly-dependent  parameters  of  the  three-point  method 
results  in  insufficient  frequency  reduction  compared  to  the  Fourier 
method.  A suitable  redefinition  of  these  parameters  gives  a three- 
point  method  that  is  stable  with  the  time  step  of  the  Fourier  methoi 
although  to  prevent  the  modified  frequency  from  becoming  zero  in  high 
latitudes  requires  excessive  computation.  Consequently  the  three- 
point  and  Fourier  methods  are  combined  to  form  a hybrid  method,  which 
is  twice  as  fast  as  the  original  Fourier  method.  An  additional 
doubling  of  the  speed  of  the  hybrid  method  is  achieved  by  a simple 
improvement  of  the  Fourier  Transform  used  in  the  Fourier  method.  The 
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hybrid  method  has  allowed  the  time  step  of  the  Rand  atmospheric  GCM 
to  be  increased  from  six  to  ten  minutes  with  a resultant  decrease  in 
time  requirement  of  37%. 
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1.  Introduction 

In  the  numerical  integration  of  the  primitive  equations  using 
explicit  time  finite  difference  schemes,  the  maximum  stable  time  step 
is  governed  by  the  highest  frequency  waves.  The  highest  frequency 
waves  that  are  admissible  in  the  primitive  system  of  equations,  the 
external  inertia-gravity  and  Lamb  waves,  have  frequencies  that 
monotonically  increase  with  decreasing  scale.  In  spherical  coordinates 
the  scale  of  any  longitudinal  wave  mode  decreases  in  the  direction  of 
the  poles  as  the  meridians  converge.  Consequently,  the  maximum  time 
step  that  is  stable  over  the  entire  sphere  (excluding  the  singular 
poles)  is  governed  by  the  smallest  scale,  highest  frequency  external 
inertia-gravity  and  Lamb  waves  near  the  poles. 

This  problem  may  be  illustrated  by  considering  a homogeneous, 
incompressible  fluid  with  a free  surface  on  a rotating  sphere. 

Ignoring  curvature  terms  except  for  the  coriolis  force,  the  governing 
equations  linearized  about  rest  are 
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where  u and  v are  respectively  the  perturbation  velocity  components  in 
the  directions  of  the  longitudinal  coordinate  X and  latitudinal  co- 
ordinate <f>,  H and  <f  are  respectively  the  mean  height  and  perturbation 
geopotential  of  the  free  surface,  g is  the  acceleration  of  gravity. 


-4- 


£1  and  a are  respectively  the  rotation  speed  and  radius  of  the  sphere, 
and  the  factor  S(n)  is  to  be  regarded  as  unity  for  the  moment.  Assuming 
normal  modes  in  the  longitudinal  direction  and  time, 


AAA 


[u,v,<J>]  (A,4>,t)  = Re{[u(<J>),v(<t>),<J>(<}>)]  exp[i (nX-wt)]} , 


where  n is  the  longitudinal  wavenumber  and  w is  the  frequency,  elimi- 

A A* 

nating  u and  4>,  and  making  the  beta-plane  approximation  about  a fixed 
latitude  4>0  (ignoring  the  variation  of  sin<}>  and  cos<j>  except  for  the 
derivative  of  the  coriolis  term)  gives 
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where 


^ ^ [„2  - -225^,  (6, 

C ' 0 ' 

and  c = (gH)^  is  the  phase  speed  of  a free  surface  gravity  wave* 
Rearranging  (6)  and  dropping  the  subscript  on  the  beta-plane  latitude 
gives  the  dispersion  relation 

1 . SSSifll  (|)2  . (2Bsin<,)2  * (|f  * l2].  (7) 


In  the  high  frequency  limit  we  obtain  the  frequency  of  the  eastward 
and  westward  propagating  inertia-gravity  waves. 
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while  in  the  low  frequency  limit  we  obtain  the  frequency  of  the  westward 
propagating  Rossby  waves, 

2 

2flnS(n)  (|) 

<2asint)2  ♦ (ff  ♦ Z2] 

As  the  latitude  of  the  beta-plane  is  moved  toward  the  poles,  cos<j>  -*■  0 

and,  for  any  wavenumber  n,  the  Rossby  wave  frequency  approaches  zero 

but  the  frequency  of  the  inertia-gravity  waves  approaches  infinity. 

Since  the  maximum  stable  time  step  for  the  explicit  integration  of  the 

primitive  equation's  is  (At)m_  = e/w,  where  e depends  on  the  particular 

max 

time  finite  difference  scheme  but  is  non-zero,  (At)max  approaches 
zero  in  the  direction  of  the  poles. 

The  dispersion  relation  (7)  is  shown  schematically  in  Fig.  1. 

For  a prescribed  time  step  At,  the  maximum  frequency  that  will  be 
stable  for  an  explicit  time  integration  scheme  is  Mmax  = e/At. 

Unless  At  is  chosen  sufficiently  small,  there  will  be  some  longitudinal 
wavenumber  nmax(<f>)  beyond  which  the  frequency  of  the  inertia-gravity 
waves  will  exceed  the  maximum  stable  frequency  |w|  , 0 Several  methods 
are  currently  employed  by  various  primitive  equation  general  circu- 
lation models  (GCM's)  to  overcome  the  need  to  use  such  a small  time 
step.  The  Geophysical  Fluid  Dynamics  Laboratory  (GFDL)  model 
(Manabe  et  al. , 1975)  applies  a Fourier  filter  every  time  step  to  each 
prognostic  variable  to  eliminate  wavenumbers  n > ^ This  method 

removes  the  small  scale,  low  frequency  Rossby  waves  along  with  the 
small  scale,  high  frequency  inertia-gravity  waves.  The  National 
Center  for  Atmospheric  Research  (NCAR)  model  (Oliger  et  al.  , 1970) 
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Wavenumber  region  where 
waves  are  either  filtered 
as  in  GFDL  model  or  are 
not  resolved  as  in  NCAR 
model  or  their  frequencies 
are  modified  as  in  UCLA 
and  Rand  models 
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Schematic  representation  of  the  dispersion  relation  (7). 
The  shading  shows  the  longitudinal  wavenumber  region 
n > n ( cfe ) where  the  frequency  of  the  inertia-gravity 

• max 

waves  exceeds  the  maximum  stable  frequency  I (tl I max  ^or  n 

prescribed  time  step  At.  The  dashed  and  dashed— dotted 
lines  illustrate  the  modified  frequencies  of  the  inertia 
gravity  and  Rossby  waves. 


decreases  the  number  of  grid  points  along  a latitude  circle  in  the 

poleward  direction  and  thus  does  not  resolve  wavenumbers  n > n (<b). 

max  Y 

The  UCLA  (Arakawa  and  Mintz,  1974)  and  Rand  (Gates  et  al. , 1971)  models 

employ  a frequency  modification  scheme  developed  by  Arakawa  for 

wavenumbers  n > nm  ($)  which  are  resolved  and  are  not  filtered.  The 

longitudinal  pressure  and  mass  flux  gradients  are  multiplied  by  a 

wavenumber-dependent  frequency  modification  factor  S(n)  as  shown  in 

(1)  and  (3)  for  the  illustrative  fluid  system.  This  modifies,  in  an 

energetically  consistent  manner,  the  influence  of  the  gradients  on  the 

tendencies;  it  does  not  modify  the  actual  fields  of  the  prognostic 

variables,  S(n)  is  defined  such  that,  for  wavenumbers  n > n (d>), 

max 

the  frequencies  given  by  (8)  for  the  inertia-gravity  waves  are  less 
than  or  equal  to  Mmax  (indicated  by  the  dashed  lines  in  Fig.  1); 
this  scheme  also  modifies  the  frequencies  of  the  Rossby  waves  for 
wavenumbers  n > nm  (<j>)  (indicated  by  the  dash-dotted  line  in  Fig.  1). 

In  this  paper  we  only  consider  models  that  employ  horizontal 
finite  differences  in  their  solution  of  the  primitive  equations.  In 
particular  we  will  focus  attention  on  grid  schemes  B,  C and  E of  the 
five  schemes  investigated  by  Winninghoff  (1968)  and  Arakawa  (1972)  for 
distributing  the  dependent  variables  on  a regular  grid  as  shown  in 
Fig.  2.  Scheme  B is  used  in  the  Rand  atmospheric  GCM  (Gates  et  al. , 
1971),  the  Goddard  Institute  for  Space  Studies  (GISS)  atmospheric  GCM 
(Somerville  et  al. , 1974)  and  the  UCLA  oceanic  GCM  (Mintz  and  Arakawa, 
1974).  Scheme  C is  used  in  the  Rand  oceanic  GCM  (Kim,  1976)  and  the 
UCLA  atmospheric  GCM  (Arakawa  and  Mintz,  1974).  Scheme  E is  used  in 
a version  of  the  Fleet  Numerical  Weather  Central  (FNWC)  atmospheric 


-9- 


GCM  (Mihok,  1974)„  In  grid  schemes  A,  B,  D and  E,  contrary  to  the 
differential  case,  the  frequency  of  the  inertia-gravity  waves  at  any 
latitude  in  two-dimensional  flow  does  not  monotonically  increase  with 
decreasing  wavelength  (Arakawa  and  Mintz,  1974).  Furthermore,  for  a 
fixed  horizontal  resolution,  the  maximum  frequency  of  the  inertia- 
gravity  waves  at  any  latitude  is  different  for  each  of  the  five  grid 
schemes.  This  means  that  the  maximum  stable  time  step  for  each 
latitude  is  also  different  for  each  grid  scheme.  However,  analogous 
to  the  differential  case,  the  frequency  of  the  inertia-gravity  waves 
grows  without  bound  as  the  poles  are  approached. 

Arakawa  has  developed  two  methods,  a "three-point"  method  (see 
Gates  et  al.3  1971)  and  a "Fourier"  method  (Arakawa,  1972),  to  modify 
the  frequencies  of  the  inertia-gravity  waves  poleward  of  a prescribed 
latitude  such  that  the  maximum  frequency  at  the  prescribed  latitude 
is  not  exceeded.  Thus  the  maximum  time  step  that  is  stable  over  the 
entire  sphere  (excluding  the  singular  poles)  is  governed  by  the 
highest  frequency  inertia-gravity  wave  at  the  prescribed  latitude. 

For  a fixed  horizontal  resolution,  this  maximum  stable  time  step  will 
be  different  for  each  of  the  five  grid  schemes.  Early  experiments 
by  Arakawa  showed  that  the  maximum  stable  time  step  for  the  Mintz- 

i 

Arakawa  two-level  GCM  was  six  minutes  with  the  three-point  method  . 
Experiments  by  the  author  showed  that  the  maximum  stable  time  step 

i 

Arakawa,  personal  communication. 
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for  this  GCM  with  the  Fourier  method  was  increased  to  ten  minutes, 
but  that  it  required  82%  more  time  than  the  GCM  with  the  three-point 
method  when  both  were  run  with  a six  minute  time  step.  Thus  the 
three-point  method  is  substantially  faster  than  the  Fourier  method 
but  the  maximum  stable  time  step  is  larger  with  the  Fourier  method 
than  with  the  three-point  method. 

The  purpose  of  the  present  study  is  to  develop  a frequency 
modification  method  that  has  the  speed  of  the  three-point  method  and 
the  maximum  stable  time  step  of  the  Fourier  method.  In  the  following, 
the  Fourier  method  is  reviewed  in  section  2.  In  section  3 the 
failure  of  the  three-point  method  to  allow  the  maximum  stable  time 
step  of  the  Fourier  method  is  analyzed.  In  section  4 it  is  shown  that 
this  failure  can  be  overcome  by  an  alternate  definition  of  the  para- 
meters of  the  three-point  method.  In  section  5 a hybrid  method  that 
combines  the  new  three-point  method  with  the  Fourier  method  is 
presented.  Finally  the  application  of  the  hybrid  method  to  several 
GCM's  is  presented  in  section  6. 
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2 Fourier  method 

The  Fourier  method  may  be  illustrated  by  considering  the  one- 
dimensional, non-rotating  case  of  (1)  - (3)  Introducing  the  coordinate 
x = a\cos4>  and  using  the  simplest  second-order  space  finite  difference 

schemes  (Arakawa,  1972),  the  one-dimensional  analogs  of  (1)  and  (3) 

« 

2 

for  grid  schemes  B,  C and  E are 

||  - - <«„*>  S(n),  ('0) 

||  = - c2  (6xu)  S(n).  (11) 


where 


(6 


a) ; 
x 'j 


Ax 


02) 


and  Ax  = aAAcos<f)  is  the  variable  longitudinal  grid  length,  AA  is  the 
angular  longitudinal  increment,  and  j is  the  grid  index  in  the 
longitudinal  direction*  Assuming  normal  modes  in  the  longitudinal 
direction  and  time. 


j 1,2,  »•.,  N 

[u,*](j,t)  = Re{[u,*]  exp[i(^P  j - cot)]},  or  (13) 

j = 1/2.3/2  ...,  N-l/2 


where  N = 2tt/AA  is  the  number  of  grid  points  along  a latitude  circle, 

A A 

and  eliminating  u or  $ gives  for  the  magnitude  of  the  frequency  of  the 


2 

Ax  and  N for  scheme  E must  be  different  from  the  corresponding 
quantities  for  schemes  A-D  if  the  total  number  of  grid  points  in  a 
given  two-dimensional  domain  is  to  be  constant. 
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eastward  and  westward  propagating  gravity  waves 

| oj | = ||sin(Trn/N)  S(n),  n = 1,  ....  N/2.  (14) 

The  frequency  modification  factor  of  the  Fourier  method  Sp(n)  is 
defined  as 


¥n>  ~ Hin['-  sin^n/N)]’  N/2’  <’5> 

where  d is  a prescribed  constant  length.,  At  (low)  latitudes  where 
Ax  > d,  Sp(n)  = 1 for  all  wavenumbers  and  the  frequency  is  not 
modified.  The  maximum  frequency 


■“Lax  = 2 c/d’ 


(16) 


occurs  at  the  latitude  where  Ax  = d.  At  (high)  latitudes  where 
Ax  < d,  Sp(n)  < 1 and  the  frequency  is  reduced  to  |oj|  for  wave- 
numbers  n > nmax(4>)  = ( N/tt ) sin  \Ax/d).  The  maximum  frequency  given 
by  (16)  is  independent  of  the  variable  longitudinal  grid  length  Ax. 

In  this  illustrative  one-dimensional  case,  the  maximum  stable  time 
step  is 


(At) 


max 


(17) 


In  the  analogous  two  dimensional  case,  the  maximum  stable  time 
step  is  approximately  governed  by  the  analog  of  (8)  for  each  grid 
scheme  at  the  latitude  where  Ax  = d„  As  illustrated  by  (17),  the 
maximum  stable  time  step  depends  on  the  choice  of  d„  At  one  extreme, 
d could  be  chosen  very  small  so  that  the  frequencies  are  modified 


d 
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only  very  near  the  poles.  At  the  other  extreme,  d could  be  chosen 
larger  than  the  maximum  Ax  so  that  the  frequencies  are  modified  at 
all  latitudes.  As  a compromise,  and  because  it  is  a natural  length 
scale  of  the  grid,  Arakawa  has  chosen  d equal  to  the  constant  latitudi- 
nal grid  length  Ay  = a A<J>  where  Aq>  is  the  angular  latitudinal  increment. 

The  Fourier  frequency  modification  factor  Sp(n)  is  shown  in  Fig.  3 
for  the  Rand  and  UCLA  atmospheric  GCM's.  Equatorward  of  38  degrees, 
Ax/A y > 1 and  Sp(n)  = 1 for  all  wavenumbers.  At  38  degrees, 

Ax/Ay  = 0.985  and  Sp(n)  < 1 for  wavenumbers  33  to  36.  As  the  poles 
are  approached,  Ax/Ay  decreases  and  Sp(n)  < 1 for  progressively  smaller 
wavenumbers.  Sp(n)  has  a minimum  value  of  0.087  (91.3%  frequency 
reduction)  at  wavenumber  36  and  86  degrees  latitude,  the  maximum 
latitude  where  the  longitudinal  pressure  gradient  and  mass  flux  are 
carried  in  these  models. 

3.  Three-point  method 

For  latitudes  where  Ax/ Ay  < 1,  the  three-point  method  is  defined 
by 


ff  = fT1 

J J 


Yf 


K-l 

j-1 


2f 


K-l 


fK-n 

J+l/’ 


J ~ 1 1 •••»  N j 

K “ If  • • • t K y 


(18) 


where  represents  the  modified  finite  difference  expression  for  the 

J 

gradient  a^/ax  or  9u/3x  at  grid  point  j at  the  K-th  iteration,  and 
y and  K are  latitude  dependent  parameters.  By  induction, 


j 


(19) 
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where  6 denotes  the  three-point  analog  to  the  second  derivative  in 
the  longitudinal  direction.  The  latitude  dependent  parameters  K and 
Y were  defined  by  Arakawa  (see  Gates  et  al.,  1971)  as 

K = [Ay/Ax],  (20) 

and 


where  [z]  denotes  the  largest  integer  less  than  or  equal  to  z„ 

We  now  analyze  the  failure  of  the  three-point  method  to  be  stable 

for  the  maximum  stable  time  step  of  the  Fourier  method,  (At) 

' 'max 

Assuming  normal  modes  in  space, 

fj  = Re{exp(i27rnj/N)},  j = 1,  ....  N,  (22) 

(19)  becomes 

f^  = [1  - 4YSin2(irn/N)]K  f.,  j = 1,  ....  N.  (23) 

J J 

The  frequency  modification  factor  of  the  three-point  method  ( n ) is 

then 

S3(n)  = [1  - 4YSin2(7rn/N)]K,  n = 1,  N/2,  (24) 

and  is  shown  in  Fig.  4 for  the  Rand  and  UCLA  atmospheric  GCM's.  Pole- 
ward  of  38  degrees  latitude,  S^(n)  < 1 for  all  wavenumbers.  S^(n) 
has  a minimum  value  of  0.0008  (99  92%  frequency  reduction)  at  86 
degrees  latitude  and  wavenumber  36.  For  the  three-point  method  to  be 
stable  for  (At)  v requires  that  S,(n)/Sc(n)  < 1 for  all  wavenumbers 

maX  J i 
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and  latitudes.  Fig,  5 shows  S3(n)/Sp(n)  for  the  Rand  and  UCLA  atmo- 
spheric GCM's,  A band  sloping  toward  decreasing  wavenumbers  with 
increasing  latitude  is  evident  where  S3(n)/Sp(n)  > 1.  In  this  band 
the  three-point  method  insufficiently  reduces  the  frequency  of  the 
inertia-gravity  waves  for  (At )max  to  be  a stable  time  step. 

On  the  basis  of  this  analysis  it  appears  that  the  failure  of  the 
three-point  method  to  be  stable  for  (At)  might  be  due  to  the 

MlaX 

definition  of  the  parameters  K and  y,  In  the  following  section, 
alternate  methods  for  determining  these  parameters  are  presented. 


4.  Alternate  parameters  for  the  three-point  method 

As  shown  above  the  frequency  modification  factors  of  the  three-point 
and  Fourier  methods  are 

S3(n)  = [1  - 4ysin2(Tm/N)]K,  n = 1,  ....  N/2,  (25) 

and 

SF(n)  = f(n) , n - 1,  ....  N/2,  (26) 


where 

f(n)  = Min  [1 ,r/sin(im/N)] , n = 1,  N/2,  (27) 

and 

r = Ax/A y - AA  cos4»/A (28) 
Defining  a function  r(n)  such  that 

Sp(n)  = [1  - 4r(n)sin2(7rn/N)]K,  n = 1 N/2,  (29) 


4 degrees  and  AA  = 5 degrees).  The 
S.(n)/S_(n)  > 1 is  shaded. 
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then  by  (26) 

. , fl/K,  . 

r(n)  = r-(n)  = ■■■*■  Ty  -An)  , n = l N/2,  (30) 

4 sin  (irn/N) 

1 /K 

where  the  positive  root  of  f is  assumed.  Since  Sp(n)  is  positive 
for  all  n,  r+(n)  is  the  only  solution  for  K odd,  while  both  r±( n)  are 
solutions  for  K even. 

The  three-point  method  will  be  stable  with  time  step  (At)  if 

max 

S3(n)/SF(n)  <1,  n = 1,  ....  N/2.  (31) 

Furthermore,  for  the  sign  of  the  modified  gradients  not  to  be  reversed 
requires  that 

S3(n)  s.  0,  n = 1 , ...»  N/2.  (32) 

By  (25),  (29)  and  (30),  conditions  (31)  and  (32)  require  for  K odd, 

T+(n)  s y £ y- . n=  1,  ....  N/2,  (33) 

4 sin  (irn/N) 

and  for  K even, 

T+(n)  s y s r"(n),  n = 1,  ....  N/2.  (34) 

Conditions  (33)  and  (34)  govern  the  range  of  y and  the  values  of  K as 
a function  of  r.  In  the  following  we  investigate  these  solutions 
separately  for  odd  and  even  K. 
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a.  K odd 

Eq0  (33)  for  K odd  may  be  written  as 

Max  r+(n)  < y < 1/4.  (35) 

The  range  of  y for  K odd  is  shown  schematically  in  Fig.  6a.  Treating 
n as  a continuous  rather  than  discrete  variable,  it  may  be  shown  that 

!1  - r1/K 

L~r—  • r a ru(K) 

(36) 

[r(K)/r]2 

4(2K+1 ) ’ r < ru^K)’ 

where 

'•„«'>  ■ (#rf-  <37> 

The  maximum  occurs  at  n = N/2  for  r > ( K)  and  at  n = 

(N/Tr)sin_1  [r/ru(K)]  for  r < ^(K).  ru(K)  is  a very  weak  monotonic 
function  for  K;  it  decreases  from  2/3  at  K = 1 to  e 1 = 0.607  as 
K -+  Substituting  (36)  into  (35)  gives 


1 - r1/K 
' 

[ru(K)/r]2] 
4(2K+1 ) ' 


s y s 1/4 


, r * ru(K) 

» r < ru(K). 


(38a) 


(38b) 


As  illustrated  in  Fig.  7,  (38a)  is  the  appropriate  solution  for 
r 2 ru(K)  while  (38b)  is  the  solution  for  rm(K)  < r < ry(K)  where 


i 


V 

0 


Eg.  (42b)  for  K even 
, Eq.(38b)for  K odd) 

r *3 


r*(K) 


4- 


. Eq  .(38a)  for  K odd 
Eq7(42a)  for  K even 


rm(K) 


V 


Range  of  ry(K) 


Fie.  7.  The  ranges  of  r where  F.qs.  (38)  are  valid 
for  K odd  and  where  Eqs.  (42)  are  valid 
for  K even.  See  text  for  the  definitions 

of  r,(K),  r (K)  and  r (K). 
i m u 


\ 


t 
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i 

: 


rjK)  = ru(K)/(2K-*-l  . 


(39) 


The  values  of  ru(K)  and  rm(K)  are  shown  in  Table  1.  Since  ru(K)  is 
bounded  by  e"'5  and  2/3,  only  (38a)  can  be  satisfied  for  r > 2/3,  only 
(38b)  can  be  satisfied  for  r < e"'5,  and  either  (38a)  or  (38b)  can  be 
satisfied  for  e”'5  < r < 2/3  depending  on  the  value  of  K.  For 
r > 2/3,  all  odd  values  of  K satisfy  (38a).  For  r < e , the  minimum 
value  of  K that  satisfies  (38b)  is  determined  by  rm(K).  For 
e”*5  ^ r < 2/3,  (38b)  is  satisfied  by  values  of  K smaller  than  some 
critical  value  while  (38a)  is  satisfied  by  K values  larger  than  this 
critical  value. 
b„  K even 

Eq.  (34)  for  K even  may  be  written  as 


Max  T+(n)  < y < Min  r”(n).  (40) 

The  range  of  y for  K even  is  shown  schematically  in  Fig.  6b.  By  (30) 
and  (27), 

1 + r1/K 

Min  T"(n)  = . (41) 

Thus  (40)  may  be  written  as 


1 - r 


1/K 


< y < 


t l/K  1 r!rU(K> 
+ r 


1 + r 


[ru(K)/r]‘ 
4(2K+1 ) 


, r < ru(K), 


(42a) 


(42b) 
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where  (36)  has  been  employed.  Eqs„  (42)  for  K even  differ  from  (38) 
for  K odd  only  in  their  right  hand  sides.  As  illustrated  in  Fig.  7, 
(42a)  is  the  appropriate  solution  for  r > ru(K)  while  (42b)  is  the 
solution  for  ^(K)  s r < ^(K)  where  r^(K)  satisfies 

rz(K)[l  + rz1/K(K )]*»  = ru(K)/(2K+l )\  (43) 

The  values  of  r (K)  and  r.(K)  are  shown  in  Table  1.  Since  r (K)  is 
u l u 

bounded  by  e 2 and  0.640  (for  K even),  only  (42a)  can  be  satisfied  for 
r > 0.640,  only  (42b)  can  be  satisfied  for  r < e , and  either  (42a) 
or  (42b)  can  be  satisfied  for  e 2 < r < 0.640  depending  on  the  value 
of  K.  For  r > 0 640,  all  values  of  K satisfy  (42a).  For  r < e"5, 
the  minimum  value  of  K that  satisfies  (42b)  is  determined  by  r^K). 

For  e 2 < r < 0.640,  (42b)  is  satisfied  by  values  of  K smaller  than 
some  critical  value  while  (42a)  is  satisfied  by  K values  larger  than 
this  critical  value. 

From  the  preceding  analysis  it  is  evident  that  for  any  r < 1, 
there  are  an  infinite  number  of  odd  K values  that  satisfy  (38)  or 
even  K values  that  satisfy  (42)  and,  for  each  value  of  K,  there  is  a 
range  of  y for  which  the  three-point  method  will  be  stable  for 
(At)max*  Clearly  additional  requirements  must  be  imposed  to  define 
a unique  solution  for  y and  K.  There  are  three  such  additional 
requirements  that  we  would  like  to  fulfill:  (a)  similarity  between 

the  dispersion  characteristics  of  the  modified  and  unmodified  inertia- 
gravity  waves,  (b)  reduction  of  the  frequency  of  the  waves  by  the 
minimum  amount  required  to  make  (At) 

max 


a stable  time  step,  and  (c) 
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selection  of  the  minimum  value  of  K so  that  the  three-point  method 
is  as  fast  as  possible  in  a numerical  solution  of  given  resolution- 
The  difficulty  of  achieving  requirement  (a)  may  be  illustrated 
by  considering  the  dispersion  relation  for  one-dimensional  gravity 
waves  shown  in  Fig„  8.  In  the  differential  case,  the  frequency 
linearly  increases  with  wavenumber;  hence  the  group  velocity  is 
constant-  In  the  finite  difference  case  without  frequency  modifi- 
cation, the  frequency  increases  with  wavenumber  more  slowly  than  the 
linear  relation  and  has  an  extremum  at  wavenumber  N/2;  hence  the 
group  velocity  monotonical ly  decreases  to  zero  at  the  highest  resolved 

wavenumber  For  (At)  to  be  a stable  time  step  when  wavenumbers 

max 

n > n are  resolved  and  are  not  filtered,  the  frequency  must  be 
max 

modified  to  be  less  than  or  equal  to  Mmax*  1°  the  Fourier  method, 
the  frequency  of  wavenumbers  n > n , is  made  equal  to  |u>|  ; hence 

their  group  velocity  is  zero  (in  two-dimensional  motion  the  longi- 
tudinal component  of  the  group  velocity  is  zero  for  these  wavenumbers). 
In  the  original  three-point  method,  the  modified  frequency  is  not  only 

too  large  to  be  stable  with  (At)  , it  exhibits  an  extremum  value 

max 

across  which  the  direction  of  the  group  velocity  becomes  opposite  to 
that  of  the  differential  and  unmodified  finite  difference  cases. 

Since  the  dispersion  relation  of  the  finite  difference  case 
without  frequency  modification  is  already  different  from  that  of  the 
differential  case  at  high  wavenumbers,  requirement  (a)  may  be  defined 
to  mean  that  the  modified  frequency  given  by  (14)  and  (25)  as 

|u>|  3 ~ sin(Tin/N)  [1  - 4ys i n2 ( nn/N ) ] K , 


(44) 
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has  no  extremum  other  than  that  at  wavenumber  N/2  due  to  the  finite 
differencing.  Again  treating  n as  a continuous  variable,  |u>|  has 
an  additional  extremum  at  wavenumber  ng  given  by 

sin  (irne/N)  - ^(2K+1 ) ° (45) 

Thus  for  |ai|  not  to  possess  an  extremum  in  1 s n < N/2  requires  that 

Y 5 4(2K+1 ) * (46) 

Imposing  this  constraint  on  (38)  for  K odd  and  on  (42)  for  K even 
gives,  for  both  K odd  and  even, 

» r ^ ru(K)  (47a) 

* Y * 4(2K+1) 

, r < ru(K).  (47b) 

Clearly  (47b)  cannot  be  satisfied,  hence  a monotonically  increasing 
frequency  is  possible  only  for  r > ru(K)  2 0,607,  For  this  range  of  r, 
requirement  (c)  can  be  fulfilled  by  choosing  the  minimum  possible  value 
of  K based  on  the  ru(K)  column  of  Table  1;  requirement  (b)  can  then  be 
fulfilled  by  choosing  y equal  to  the  left  side  of  (47a),  For  r < 0,607, 
either  the  Fourier  method  must  be  employed  to  avoid  the  extremum  or 
requirement  (a)  must  be  redefined. 

Requirement  (a)  may  be  defined  to  mean  that  the  modified  frequency 
given  by  (44)  may  possess  an  additional  extremum  but  may  not  equal 
zero  at  any  wavenumber  (this  definition  was  implicit  in  the  original 


1 - r1^  ) 
* 

[ru(K)/r]2J 
4(2K+1)  ' 
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three-point  method)  By  (44),  |w|  equals  zero  at  wavenumber  nQ  given 

by 


sin2(TTno/N)  = ^ 


(48) 


For  M not  to  equal  zero  in  1 < n < N/2  thus  requires  that 


Y < 1/4.. 


(49) 


Imposing  this  constraint  on  (38)  for  K odd  and  on  (42)  for  K even 
gives  for  both  K odd  and  even. 


4 


[ru(K)/r]2 

4(2K+1) 


< Y < 1/4 


, r > ru(K) 


. r < ru(K), 


(50a) 


(50b) 


Requirement  (c)  may  be  fulfilled  by  choosing  K = 1 for  r > 0.667  and 
the  minimum  possible  K under  the  rm(K)  column  of  Table  1 for  r < 0.667. 
Note  that  for  small  r,  by  (39),  K a r ■,  hence  as  the  poles  are 
approached  many  iterations  are  required  to  prevent  the  frequency  from 
equaling  zero  (for  discrete  n the  frequency  will  not  equal  zero  if 
nQ  * integer  but  will  be  very  small  for  wavenumbers  near  nQ),  Require- 
ment (b)  may  then  be  satisfied  by  choosing  y equal  to  the  left  hand 
side  of  (50a)  and  (50b)  for  r > 0.667  and  r < 0.667,  respectively. 

If  requirement  (a)  is  dropped  altogether,  the  fastest  version  of 
the  three-point  method  can  be  obtained  from  Table  1 as  follows.  From 
the  r(j(K)  column  we  see  that  K = 1 is  sufficient  with  (38a)  for 
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r > 0.667.  For  this  range  of  r,  the  modified  frequency  will  not  equal 

zero  and  will  be  monotonically  increasing.  Requirement  (b)  may  then  be 

satisfied  by  choosing  y equal  to  the  left  hand  side  of  (38a).  From  the 

r (K)  column  we  see  that  K = 1 is  sufficient  with  (38b)  for  0 385 
m 

<,  r < 0.667.  For  this  range  of  r,  the  modified  frequency  will  not 
equal  zero.  Requirement  (b)  may  then  be  satisfied  by  choosing  y equal 
to  the  left  hand  side  of  (38b).  By  comparing  the  r^(K)  and  rm(K) 
columns  we  see  that  r^(K)  < rm(K+l)  where  K is  even.  Consequently, 
for  r < 0.385,  we  may  choose  the  minimum  even  K based  on  the  r^(K) 
column  and  select  y from  (42b);  the  precise  choice  of  y determines 
where  the  frequency  achieves  an  extremum  (by  45))  and  where  it  equals 
zero  (by  (48)). 

Since  it  is  undesirable  for  the  modified  frequency  to  equal  zero 
at  any  wavenumber,  and  because  the  number  of  iterations  K required  to 
prevent  this  rapidly  increases  as  the  poles  are  approached,  the  hybrid 
method  presented  in  the  following  section  has  been  developed., 


5.  Hybrid  method 

In  this  section  the  three-point  and  Fourier  methods  are  combined 
to  form  a hybrid  method.  The  frequency  modification  factor  of  the 
hybrid  method  S^(n)  is  defined  as 


SH(n) 


S3(n), 

SF(n), 


r > rmin^max^’ 


r < rmin^max^  * 


(51) 


where  S3(n)  is  given  by  (25),  Sp(n)  is  given  by  (26)  and  (27),  and 
the  meaning  of  rm-j n ( Kmax ) 1S  explained  below. 

As  shown  in  section  4,  the  right  side  of  (50)  is  the  necessary 
condition  for  S3(n)  * 0.  Since  y < 1/4,  S3 ( n ) is  a monotonically 
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decreasing  function  of  y„  To  fulfill  requirement  (b)  we  may  therefore 
choose  the  minimum  value  of  y given  by  (50)  as 


r * ru(K) 


r(K) 

m 


< r < ry(K) , 


(52) 


where  (39)  has  been  used*  Since  (52)  is  identical  to  (35), 

Y = Max  T+(n)  = r+(m)  where  m = N/2  for  r > ry(K)  and  m = 

(N/tt)  sin_1[r/ru(K)]  for  r < ru(K)„  By  the  definition  of  r+(n), 

S3(n)/SF(n)  =1  for  n = m and  S3(n)/Sp(n)  <1  for  n / m. 

For  a given  value  of  r,  (52)  prevents  S3(n)  = 0 by  requiring  that 

rm(K)  < r,  By  (39)  this  requires  that  K a r"2  + « as  r -*■  0 towards 

the  poles.  At  some  value  of  r it  therefore  becomes  more  economical 

to  switch  from  the  three-point  method  to  the  Fourier  method.  However, 

while  S3 ( n ) = 0 is  prevented  by  (52),  the  minimum  value  of  S3(n)/Sp(n) 

over  n = 1,  ....  N/2  becomes  quite  small  as  r -►  r (K)  for  fixed  K. 

m 

Letting  M represent  this  minimum  value,  it  may  be  shown  that 


M = 


1 - (1  - r1/K) 

rf. 

r ^ ry(K) 

(53a) 

1 - rjm 

K 

rp(K)  s r < ru(K) 

(53b) 

1 - (rm(K)/r)2 

s 

r * 

rm(K)  < r s rp(K), 

(53c) 

* 
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where 

rp<K>  ' rro(K)j'  - -p,/K(K)[.  - rm2(K)][  A (54) 

M is  shown  in  Fig.  9 as  a function  of  r for  several  values  of  K.  For 
any  K,  M decreases  from  1 at  r = 1 to  [1  -r  2 ( K) ] K at  r = r (K),  is 
constant  for  rp(K)  < r < ^(K),  and  then  approaches  zero  as  r ■*  rm(K). 
Therefore,  to  satisfy  requirement  (b),  we  impose  a lower  bound  on  M 

M > 6,  r < r (K),  (55) 

where  0 < 6 < e^e  = .832  (the  upper  bound  on  6 is  given  by  (53b)  as 
K ■*  °°).  As  shown  in  Fig.  9,  imposing  a nonzero  6 increases  the 
minimum  value  of  r for  a given  K.  Letting  rmin(K)  denote  this  minimum 
value,  (55)  and  (53c)  give 

' ''nM1  - Ki,,<K>)'/Kr.  <56) 

Fig.  10  shows  rm^n(K)  for  several  values  of  6.  For  nonzero  6, 

rmi n ( 1^)  decreases  substantially  with  increasing  K only  to  about  K = 10 

and  thereafter  is  quasi-constant.  Consequently,  for  a given  horizontal 

grid  resolution  and  nonzero  5,  there  is  an  effective  maximum  value  of 

K,  K , beyond  which  no  higher  latitudes  may  be  treated  by  the  three- 

point  method.  For  example,  for  the  Rand  and  UCLA  atmospheric  GCM's 

with  AiJ)  = 4 degrees,  AX  = 5 degrees  and  6 = 0.8,  r = .428  at  4>  = 70 

degrees  can  be  treated  by  K = 5 but  r =,,345  at  = 74  degrees  cannot 

be  treated  by  any  value  of  K.  For  r < r . (K  ) the  Fourier  method 

min  max 


must  be  used. 


min 
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Fig.  10.  The  minimum  value  of  r as  a function  of  K for 
several  values  of  6. 
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The  parameters  K and  y for  the  three-point  method  can  be  deter- 
mined from  (52)  as  a function  of  r.  To  fulfill  requirement  (c),  we 
choose  the  minimum  value  of  K that  is  consistent  with  (52)  and  (56). 
This  gives  the  following  selection  rule  for  y: 

, r 2 ru(l)  » .667  (57a) 

)[rJD/r]2 

y-  S . rm,n(l)  s r < ru(l)  (57b) 

[rJKJ/r]2 

V —4—  ’ WK>  s r < K ' 2-  W <57c> 

6.  Application  to  GCM's 

As  an  illustration  we  present  the  hybrid  method  for  the  Rand  and 

UCLA  atmospheric  GCM's  (A<j>  = 4 degrees  and  AA  = 5 degrees).  We  choose 

6 = .8  to  minimize  the  frequency  reduction  of  the  three-point  method. 

The  parameters  K and  y are  shown  in  Table  2 as  a function  of  latitude 

along  with  the  corresponding  value  of  r and  the  appropriate  member  of 

(57).  At  34  degrees  latitude,  r > 1,  hence  frequency  modification  is 

not  required.  From  38  degrees  to  54  degrees  latitude,  r > .667,  hence 

K = 1 and  y is  given  by  (57a).  From  58  degrees  to  66  degrees  latitude, 

r > r ^ (1)  = .495,  hence  K = 1 and  y is  given  by  (57b).  At  70  degrees 

latitude,  r . (5)  < r < r . (4),  hence  K = 5 and  y is  given  by  (57c). 
min  min 

At  74  degrees  latitude,  r < rml- n(°°) , hence  this  and  higher  latitudes 
must  be  treated  by  the  Fourier  method. 

The  frequency  modification  factor  of  the  hybrid  method  S^(n)  is 
shown  in  Fig.  11  for  the  Rand  and  UCLA  atmospheric  GCM's.  Comparison 
with  Fig.  4 for  Sp(n)  and  Fig.  5 for  S^(n)  shows  that  the  only 
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Table  2.  Three-point  method  parameters  K and  y for  the 

Rand  and  UCLA  atmospheric  GCM's.  A<J>  = 4 degrees, 
AX  = 5 degrees  and  6 = „8. 


Latitude 

r 

. 

Eq.  (57) 

K 

100  y 

34 

1.036 

NOT  REQUIRED 

38 

.985 

a 

1 

.375 

42 

.929 

a 

1 

1.78 

46 

.868 

a 

1 

3.30 

50 

.803 

a 

I 

4.92 

54 

.735 

a 

1 

6.64 

58 

.662 

b 

1 

8.45 

62 

.587 

b 

1 

10.8 

66 

.508 

b 

1 

14.4 

70 

.428 

c 

5 

4.80 

74 

.345 

NOT  POSSIBLE 
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qualitative  difference  is  that  Su(n)  increases  between  70  degrees 
and  74  degrees  latitude  for  wavenumbers  1 through  9.  This  occurs 
because  Sc(n)  e 1 for  n < nm  (cf>)  but  S~(n)  < 1.  Fig..  12  shows  that 
.8  < S^(n)/Sp-(n)  < 1 over  the  entire  latitude-wavenumber  domain,. 

The  time  requirements  of  several  frequency  modification  methods 
relative  to  the  original  three-point  method  are  shown  in  Table  3 for 
the  Rand  and  UCLA  GCM's.  The  original  Fourier  method  uses  a classical 
Fourier  Transform  (see  Appendix)  and  is  15.2  times  slower  than  the 
original  three-point  method  (the  Rand  atmospheric  GCM  with  this  original 
Fourier  method  requires  82 % more  time  than  the  GCM  with  the  original 
three-point  method  when  both  are  run  with  a six  minute  time  step),.  The 
hybrid  method  that  combines  the  three-point  method  with  parameters 
given  in  Table  2 and  the  original  Fourier  method  reduces  the  time 
requirement  by  a factor  of  two„  To  further  reduce  the  time  requirement, 
the  classical  Fourier  Transform  was  replaced  by  a Fast  Fourier  Trans- 
form. As  shown  in  Table  3,  the  Fourier  method  with  the  Fast  Fourier 
Transform  is  slower  than  the  Fourier  method  with  the  classical  Fourier 
Transform.  This  occurs  because  the  Fast  Fourier  Transform  transforms 
all  wavenumbers  at  a particular  latitude  while  the  classical  Fourier 
Transform  only  transforms  the  required  wavenumbers  n > nmax(4))  (see 
Appendix).  However,  the  hybrid  method  with  the  Fast  Fourier  Transform 
is  slightly  faster  than  the  hybrid  method  with  the  classical  Fourier 
Transform  because,  at  the  (high)  latitudes  where  the  Fourier  method 
is  used,  nearly  all  wavenumbers  must  be  transformed.  As  shown  in  the 
Appendix,  the  speed  of  the  classical  Fourier  Transform  can  be  nearly 
doubled  by  a simple  improvement.  Table  3 shows  that  the  Fourier  method 
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Table  3o  Approximate  time  requirements  of 
several  frequency  modification 
methods  relative  to  the  original 
three-point  method.  A<p  = 4 degrees, 
AA  = 5 degrees  and  6 = „8. 


METHOD 

TIME 

Original  Three-Point  Method 

1 

Fourier  method  with  classical  Fouri er  Transform 

15.2 

Hybrid  method  with  classical  Fourier  Transform 

7.7 

Fourier  Method  with  Fast  Fourier  Transform 

19.7 

Hybrid  Method  with  Fast  Fourier  Transform 

6,5 

Fourier  Method  with  Improved  Fourier  Transform 

8*2 

Hybrid  Method  with  Improved  Fourier  Transform 

3.8 
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with  this  improved  Fourier  Transform  is  almost  twice  as  fast  as  the 
original  Fourier  method  and  that  the  hybrid  method  with  the  improved 
Fourier  Transform  is  four  times  faster  than  the  original  Fourier  method. 
Thus,  regardless  of  the  Fourier  Transform  used,  the  hybrid  method  is 
twice  as  fast  as  the  Fourier  method. 

Finally,  Table  4 shows  the  percentage  decrease  in  GCM  running  time 
for  several  of  the  models  that  are  using  the  hybrid  method  with  the  improved 
Fourier  Transform.  The  Rand  atmospheric  GCM  previously  used  the  origi- 
nal three-point  method  with  its  maximum  stable  time  step  of  six  minutes. 

The  hybrid  method,  while  3.8  times  slower  than  the  original  three-point 
method,  allows  the  time  step  to  be  increased  to  the  ten  minute  maximum 
stable  time  step  of  the  Fourier  method.  This  results  in  a 37%  decrease 
in  the  overall  running  time.  In  the  other  models,  the  8.6  - 33%  de- 
crease in  overall  running  time  results  from  the  fourfold  increase  in 
speed  of  the  hybrid  method  compared  to  the  original  Fourier  method. 

(The  difference  in  the  percentage  decrease  in  running  time  for  these 
GCM's  is  due  to  differences  in  their  time  marching  procedures  and  in 
their  parameterizations  of  physical  processes.) 

7.  Summary 

Arakawa  has  developed  two  methods,  a "three-point"  method  and  a 
"Fourier"  method  to  modify  the  frequencies  of  inertia-gravity  waves 
poleward  of  a prescribed  latitude  such  that,  without  Fourier  filtering 
or  reducing  the  longitudinal  resolution,  these  waves  satisfy  the 
linear  stability  criterion  based  on  the  maximum  stable  time  step  at 
the  prescribed  latitude.  Experiments  have  shown  that  the  three-point 


Table  4.  Percentage  decrease  in  GCM  running  time  for  several  models 
using  the  hybrid  method  with  improved  Fourier  Transform. 
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method  is  substantially  faster  than  the  Fourier  method  but  that  the 
maximum  stable  time  step  is  larger  with  the  Fourier  method  than  with 
the  three-point  method . 

An  analysis  of  the  three-point  method  shows  that  the  modified 
frequencies  of  the  inertia-gravity  waves  are  too  large  to  be  stable 
with  the  maximum  stable  time  step  of  the  Fourier  method,  and  suggests 
that  this  insufficient  frequency  reduction  is  due  to  the  definition 
of  the  latitude-dependent  parameters  K and  y of  the  three-point  method. 
It  is  shown  that  these  parameters  can  be  redefined  in  such  a way  that 
the  three-point  method  is  stable  with  the  maximum  stable  time  step  of 
the  Fourier  method;  however,  this  stability  requirement  alone  is  not 
sufficient  to  determine  unique  values  of  K and  y as  a function  of 
latitude.  These  parameters  can  be  uniquely  determined  by  imposing 
the  additional  requirements  of:  (a)  similarity  between  the  dispersion 

characteristics  of  the  modified  and  unmodified  inertia-gravity  waves, 
(b)  reduction  of  the  frequency  of  these  waves  by  the  minimum  amount 
required  for  stability,  and  (c)  selection  of  the  minimum  number  of 
iterations  K so  that  the  three-point  method  is  as  fast  as  possible. 

It  is  shown  that  the  resultant  modified  frequencies  of  the  inertia- 
gravity  waves  do  not  monotonical ly  increase  with  wavenumber  in  high 
latitudes  and  that  y < 1/4  is  required  to  prevent  the  modified 
frequencies  from  becoming  zero. 


The  constraint  y < 1/4  requires  that  K ■+•  °° 


as  the  poles  are 
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approached,  At  some  latitude  it  therefore  becomes  more  economical  to 
switch  from  the  three-point  method  with  new  parameters  to  the  Fourier 
method.  Thus  a hybrid  method  is  defined  that  uses  the  three-point 
method  at  low  latitudes  and  the  Fourier  method  at  high  latitudes,  A 
selection  rule  is  given  for  K and  y such  that  the  modified  frequency 
of  the  three-point  method  is  not  less  than  a prescribed  fraction  of 
the  modified  frequency  that  would  be  given  by  the  slower  Fourier  method. 

The  hybrid  method  is  stable  with  the  maximum  stable  time  step  of 
the  Fourier  method  and  is  twice  as  fast.  An  additional  doubling  of 
the  speed  of  the  hybrid  method  is  achieved  by  a simple  improvement 
of  the  Fourier  Transform  used  in  the  Fourier  method. 

The  hybrid  method  has  allowed  the  time  step  of  the  Rand  atmo- 
spheric GCM  to  be  increased  from  six  to  ten  minutes  with  a resultant 
decrease  in  time  requirement  of  37%.  The  fourfold  increase  in  speed 
of  the  hybrid  method  compared  to  the  original  Fourier  method  has 
resulted  in  a decrease  in  time  requirement  of  8.6%  in  the  UCLA 
12-level  atmospheric  GCM  and  33%  in  the  Mintz-Arakawa  3-level  model 
for  Mars. 
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APPENDIX 


Improved  Fourier  Transform 


The  longitudinal  pressure  and  mass  flux  gradients  at  longitudinal 

grid  point  j may  be  represented  by  a Fourier  series 

1 N/2  N/2-1  , 

f-i  = 2 A(0)  + S A(n)  cos  (■=£■  nj ) + 2 B(n)  sin  (^  nj) , 

J n=l  " n=l  H 

j = 1,  ....  N.  (Al) 

In  the  Fourier  method,  fj  is  modified  to 

* 1 N/2 

f,  = j A( 0 ) + I SF(n)  A(n)  cos  nj ) + 

J c n=l  h " 

N/|_1  Sc(n)  B(n)  sin  (-2-£  nj)  , j = 1,  ....  N,  (A2) 
n=l  h " 

where  Sp(n)  is  given  by  (15).  Subtracting  (Al)  from  (A2)  gives 

* N/2  2TT 

f . = f . - I [1  - S,(n)]  A(n)  cos  nj)  + 

J J n=l  h 


N/2-1  ?.T 

2 [1  - S F ( n ) ] B(n)  sin  nj)],  j = 1,  ...,  N.  (A3) 

n=l 

Since  Sc(n)  = 1 for  n < n (<{>),  the  lower  bound  of  the  summations  in 
F'  ' max 

(A3)  should  be  replaced  by  nm:,  (<f>)  + 1.  In  the  original  Fourier  method, 

max 

the  coefficients  A(n)  and  B(n)  are  calculated  by  the  classical  Fourier 
Transform 


A(n) 


= Q(n) 


N 

j=l 


f.  cos  (— -nj),  n = nmax(0)  + 1,  N/2, 


B(n ) 


2 

N 


N 9TT 

2 sin  (^rr  nj), 
j-1  J N 


" - "max(*>  + 1 N/2  - 1 


(A4) 
( A5) 
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where 


a(n)  = 


2 , n * N/2 


1 , n = N/2. 


EqSo  (A4)  and  (A5)  together  require  2N  multiplications  and  2N-2 
additions  for  each  n. 

The  speed  of  this  Fourier  Transform  can  be  improved  as  follows. 
Writing  (A4)  and  (A5)  as 


IN  / c 

A(n)  = 2 fj  cos  nj)  + fj+N/2  cos  n(j  + N/2)], 

J • 


9 L 9tt  ?tt 

B(n)  = ^ ^ sin  nj)  + fj+N/2  sin  nU  + N/2)]. 


and  using  the  trigonometric  relations 


cos  n(j  + N/2)]  = (-  l)n  cos  (-2-J  nj), 

sin  n(j  + N/2)]  = (-  1)"  sin  nj). 


gives 


A(n)  * SM  gj(n)  cos  nj). 


o IV  c ? 

B(n)  = £ I g • ( n ) sin  (*J  nj), 
N j=l  J 


where 


fj  ~ f j+N/2 , n odd 


9j(n)  - 


j = 1,  N/2  (A9) 


f,i  + f j+N/2,  n even. 
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Eqs<.  (A7)  and  (A8)  together  require  N multiplications  and  N-2  additions 

for  each  n,  and  the  calculation  of  the  N different  g.(n)  requires  N 

0 

additions.  The  improved  Fourier  Transform  is  thus  almost  twice  as  fast 
as  the  classical  Fourier  Transform. 


-5507  A FAST  NUMERICAL  METHOD  FOR  EXPLICIT  INTEGRATION  OF  THE  Schlesinger 

PRIMITIVE  EQUATIONS  NEAR  THE  POLES 


